{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cudf\n",
    "import numpy as np\n",
    "from cuml.datasets import make_blobs\n",
    "from cuml.neighbors import NearestNeighbors as cuNearestNeighbors\n",
    "from sklearn.neighbors import NearestNeighbors as skNearestNeighbors\n",
    "import argparse, numpy as np, pandas as pd, timeit, os, subprocess\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "def load_LAS(las_fname, dtype='float32'):\n",
    "    \"\"\"\n",
    "    Load LAS or LAZ file (only coordinates) and return pc_xyz and xy vectors. Converts float64 to float32 by default, unless you set dtype='float64'\n",
    "    \"\"\"\n",
    "    from laspy.file import File\n",
    "\n",
    "    inFile = File(las_fname, mode='r')\n",
    "    pc_pc_xyz = np.vstack((inFile.get_x()*inFile.header.scale[0]+inFile.header.offset[0], inFile.get_y()*inFile.header.scale[1]+inFile.header.offset[1], inFile.get_z()*inFile.header.scale[2]+inFile.header.offset[2])).transpose()\n",
    "\n",
    "    #setting datatype to float32 to save memory.\n",
    "    if dtype == 'float32':\n",
    "        pc_pc_xyz = pc_pc_xyz.astype('float32')\n",
    "    return pc_pc_xyz\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading input file: test_data/Pozo_WestCanada_clg.laz... loaded 3,348,668 points\n"
     ]
    }
   ],
   "source": [
    "#Load in data\n",
    "lasfname = 'test_data/Pozo_WestCanada_clg.laz'\n",
    "print('Loading input file: %s... '%lasfname, end='', flush=True)\n",
    "pc_xyz = load_LAS(lasfname, dtype='float32')\n",
    "mean_x = np.nanmean(pc_xyz[:,0])\n",
    "mean_y = np.nanmean(pc_xyz[:,1])\n",
    "mean_z = np.nanmean(pc_xyz[:,2])\n",
    "#subtracting mean value to center around 0 - UTM coordinates tend to be large!\n",
    "pc_xyz[:,0] = pc_xyz[:,0] - mean_x\n",
    "pc_xyz[:,1] = pc_xyz[:,1] - mean_y\n",
    "pc_xyz[:,2] = pc_xyz[:,2] - mean_z\n",
    "print('loaded %s points'%\"{:,}\".format(pc_xyz.shape[0]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using cKDTree to calculate neighborhood. Leaf sizes of ~20 have been shown to work well with lidar point clouds (cf. [https://github.com/UP-RS-ESP/LidarPC-KDTree](https://github.com/UP-RS-ESP/LidarPC-KDTree)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pc_generate_cKDTree(pc_xyz, leafsizei=10):\n",
    "    try:\n",
    "        from scipy import spatial\n",
    "    except ImportError:\n",
    "        raise pc_generate_cKDTree(\"scipy not installed.\")\n",
    "    pc_xyz_cKDTree_tree = spatial.cKDTree(pc_xyz, leafsize=leafsizei)\n",
    "    return pc_xyz_cKDTree_tree\n",
    "\n",
    "def pc_query_cKDTree(pc_xyz_cKDTree_tree, pc_xyz, k=10):\n",
    "    pc_cKDTree_distance, pc_cKDTree_id = pc_xyz_cKDTree_tree.query(pc_xyz, k=k, n_jobs=-1)\n",
    "    return pc_cKDTree_distance, pc_cKDTree_id\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pc_density_sphere(d, k=10):\n",
    "    dens =  k / np.pi / d[:, -1]**2\n",
    "    return dens\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pc_density_ellipsoid(pc_xyz, pc_ids, k=10):\n",
    "    #estimate density in x, y, z directions using an ellipsoid and returning lengthening ratios in x, y, z directions\n",
    "\n",
    "    x = np.max(pc_xyz[pc_ids[:,1::],0],axis=1) - np.min(pc_xyz[pc_ids[:,1::],0],axis=1)\n",
    "    y = np.max(pc_xyz[pc_ids[:,1::],1],axis=1) - np.min(pc_xyz[pc_ids[:,1::],1],axis=1)\n",
    "    z = np.max(pc_xyz[pc_ids[:,1::],2],axis=1) - np.min(pc_xyz[pc_ids[:,1::],2],axis=1)\n",
    "    dens_ellipsoid = k / (4/3 * np.pi * x * y * z )\n",
    "    pts_x_y_ratio = x / y #values < 1 indicate X axis is compressed (and Y axis is longer)\n",
    "    pts_x_z_ratio = x / z #values < 1 indicate X axis is compressed (and Z axis is longer), values > 1 indicate X axis is longer\n",
    "    pts_y_z_ratio = y / z #values < 1 indicate X axis is compressed (and Z axis is longer), values > 1 indicate X axis is longer\n",
    "    return dens_ellipsoid, pts_x_y_ratio, pts_x_z_ratio, pts_y_z_ratio\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
